Temps = [0.35,0.35,0.45,0.6,0.75,1.3,3.75,7,9.6];
TempsErr = [0.1,0.1,0.1,0.15,0.1,0.3,0.5,0.5,0.3];
%TODO: get reasonable temp error bars

Ds = [8,6.4,3.6,2.5,3.5,2.7,2.4,5.7,8.5];
DsErr = [2.1,1,0.6,0.9,1,0.4,1.1,2.0,3.2];
Gammas = [4,4.3,6,6.1,7,5.8,7.2,10.4,11.8];
GammasErr = [0.9,0.5,0.6,0.7,0.9,0.7,1.1,3.1,3.5];
SoundSpeed = [10.7,11.6,11.8,12.2,11.2,9.5,9.6,8.7,7.5];
SoundSpeedErr = [0.5,0.3,0.3,0.4,0.4,0.2,0.5,0.9,2];
Visc = zeros(size(Temps));
ViscErr = zeros(size(Temps));

FullD = Ds + (SoundSpeed.^2./Gammas);
SoundSpdOverGamma_Err = SoundSpeed.^2./Gammas.*sqrt((2*SoundSpeedErr./SoundSpeed).^2 + (GammasErr./Gammas).^2);
FullDErr = sqrt(DsErr.^2 + SoundSpdOverGamma_Err.^2);

NRows = 2;
NCols = 3;

figure;
subplot(NRows,NCols,1)
errorbar(Temps,Ds,DsErr,DsErr,TempsErr,TempsErr);
grid on;
xlabel('T (t)');
ylabel('D_{Inc} (latt^2/ms)');
ax = gca;
ax.YLim(1) = 0;
title('Incoherent/Microscope D');

subplot(NRows,NCols,2)
errorbar(Temps,Gammas,GammasErr,GammasErr,TempsErr,TempsErr);
grid on;
xlabel('T (t)');
ylabel('\Gamma (KHz)');
ax = gca;
ax.YLim(1) = 0;
title('Momentum Relaxation Rate');

subplot(NRows,NCols,3)
errorbar(Temps,Visc,ViscErr,ViscErr,TempsErr,TempsErr);
grid on;
xlabel('T (t)');
ylabel('\nu (latt^2/ms)');
ax = gca;
ax.YLim(1) = 0;
title('Viscosity');

subplot(NRows,NCols,4)
errorbar(Temps,SoundSpeed,SoundSpeedErr,SoundSpeedErr,TempsErr,TempsErr);
grid on;
xlabel('T (t)');
ylabel('c (latt/ms)');
ax = gca;
ax.YLim(1) = 0;
title('Sound Speed');

subplot(NRows,NCols,5)
errorbar(Temps,FullD,FullDErr,FullDErr,TempsErr,TempsErr);
grid on;
xlabel('T (t)');
ylabel('D (latt^2/ms)');
ax = gca;
ax.YLim(1) = 0;
title('D = c^2/\Gamma + D_{inc}');

subplot(NRows,NCols,6)
DInv = 1./FullD;
DInvErr = FullDErr./FullD.^2;
errorbar(Temps,DInv,DInvErr,DInvErr,TempsErr,TempsErr);
hold on;
xlabel('T (t)')
ylabel('1/D ms/latt^2');
grid on;
ax = gca;
ax.YLim(1) = 0;
title('1/D');

suptitle('Model Params Vs. T')

figure;
subplot(1,2,1)
errorbar(Temps,Ds,DsErr,DsErr,TempsErr,TempsErr);
hold on;
errorbar(Temps,SoundSpeed.^2./Gammas,SoundSpdOverGamma_Err,SoundSpdOverGamma_Err,TempsErr,TempsErr);
grid on;
xlabel('T (t)')
ylabel('ms/latt^2');
grid on;
ax = gca;
ax.YLim(1) = 0;
legend({'D_{Inc}','c^2/\Gamma'});
title('Coherent-ness')

subplot(1,2,2)
Ratio = Ds./FullD;
RatioErr = Ds./FullD.*sqrt((DsErr./Ds).^2+(FullDErr./FullD).^2);
errorbar(Temps,Ratio,RatioErr,RatioErr,TempsErr,TempsErr);
grid on;
ax = gca;
ax.YLim(1) = 0;
xlabel('T (t)');
ylabel('D_{Inc}/D');
